##############
### Tables ###
##############


# Table 1
coef_table1 <- stargazer(model1, model1.re, model3, model3.re, model5, model5.re,
                         column.labels=c("log Intra-region trade", "log Inter-region trade", 
                                         "Region trade proportion"),
                         column.separate = c(2, 2, 2),
                         model.names=FALSE, style="aer",
                         ci=FALSE, p.auto=TRUE,
                         dep.var.labels.include=FALSE,
                         dep.var.labels="",
                         table.placement="htbp",
                         font.size="footnotesize",
                         model.numbers=TRUE, #float.env="sidewaystable",
                         notes.label=" ",
                         digits=2, notes="", df=FALSE, omit.stat=c("rsq", "aic", "bic", "ll"),
                         order = c(1, 2, 55, 3, 8, 4, 5, 6, 7),
                         omit = "fyear",
                         covariate.labels = c("Regional power NHHI", "Regional power Polity", "NHHI X Reg. power Polity",
                         "Region cinc total", "Count of states in region", "log GDP", "log Population", "Polity combined score", "Armed conflict"),
                         title = "Concentration of regional military capabilities and state-level trade. Note: odd-numbered models include fixed effects by state and
                         year; even numbered models include random effects (intercepts) by region and state, and fixed effects by year.",
                         label = "table:coefficients1",
                         notes.align = "r",
                         star.cutoffs = c(0.05, 0.01, 0.001))
                        
write(x=coef_table1, file=paste0(getwd(), "/tables/Coef-Table1.tex"))


# Appendix Table 1
coef_tableA1 <- stargazer(model2, model2.re, model4, model4.re, model6, model6.re,
                          column.labels=c("log Intra-region trade", "log Inter-region trade", 
                                          "Region trade proportion"),
                          column.separate = c(2, 2, 2),
                          model.names=FALSE, style="aer",
                          ci=FALSE, p.auto=TRUE,
                          dep.var.labels.include=FALSE,
                          dep.var.labels="",
                          table.placement="htbp",
                          font.size="footnotesize",
                          model.numbers=TRUE, #float.env="sidewaystable",
                          notes.label=" ",
                          digits=2, notes="", df=FALSE, omit.stat=c("rsq", "aic", "bic", "ll"),
                          order = c(1, 2, 55, 3, 8, 4, 5, 6, 7),
                          omit = "fyear",
                          covariate.labels = c("Regional power NHHI", "Regional power Polity", "NHHI X Reg. power Polity",
                                               "Region cinc total", "Count of states in region", "log GDP", "log Population", "Polity combined score", "Armed conflict"),
                          title = "Concentration of regional military capabilities and state-level trade. Note: odd-numbered models include fixed effects by state and
                          year; even numbered models include random effects (intercepts) by region and state, and fixed effects by year.",
                          label = "table:coefficientsA1",
                          notes.align = "r",
                          star.cutoffs = c(0.05, 0.01, 0.001))

write(x=coef_tableA1, file=paste0(getwd(), "/tables/Coef-TableA1.tex"))

# Appendix Table 2 - conflict-excluded models
coef_tableA2 <- stargazer(model1nc, model2nc, model3nc, model4nc, model5nc, model6nc,
                         column.labels=c("Intra-region trade", "Inter-region trade", 
                                         "Intra-region proportion"),
                         column.separate = c(2, 2, 2),
                         model.names=FALSE, style="aer",
                         ci=FALSE, p.auto=TRUE,
                         dep.var.labels.include=FALSE,
                         dep.var.labels="",
                         table.placement="t",
                         font.size="footnotesize",
                         model.numbers=TRUE, #float.env="sidewaystable",
                         notes.label=" ",
                         digits=2, notes="", df=FALSE, omit.stat=c("rsq"),
                         order = c(1, 2, 7, 3, 4, 5, 6),
                         covariate.labels = c("Regional power NHHI", "Regional power Polity", "NHHI X Reg. power Polity",
                                              "Region cinc total", "log GDP", "log Population", "Polity combined score"),
                         title = "Concentration of regional military capabilities and state-level trade, excluding variable for armed conflict. 
                         Note: odd-numbered models include fixed effects by state and year; even numbered models include random 
                         effects (intercepts) by region and state, and fixed effects by year.",
                         label = "table:coefficientsA2",
                         notes.align = "r",
                         star.cutoffs = c(0.05, 0.01, 0.001))

write(x=coef_tableA2, file=paste0(getwd(), "/tables/Coef-TableA2.tex"))

# Appendix Table 3 - dyadic models
coef_tableA3 <- stargazer(model1dr, model1dry, model2dr, modelP1, modelP2,
                          column.labels=c("Random Intercept Models", "Quasi-poisson Models"),
                          column.separate = c(3, 2),
                          model.names=FALSE, style="aer",
                          ci=FALSE, p.auto=TRUE,
                          dep.var.labels.include=FALSE,
                          dep.var.labels="",
                          table.placement="htbp",
                          font.size="footnotesize",
                          notes.label=" ",
                          model.numbers=TRUE, #float.env="sidewaystable",
                          digits=2, notes="", df=FALSE, omit.stat=c("rsq", "aic", "bic"),
                          order = c(1, 2, 15, 3, 4, 5, 6, 16, 7, 8, 9, 10, 11, 12, 13, 14),
                          covariate.labels = c("Exporter regional NHHI", "Exporter power Polity", "Exporter NHHI X Polity",
                                               "Exporter CINC total", "Exporter state count", "Importer regional NHHI", "Importer power Polity", "Importer NHHI X Polity",
                                               "Importer CINC total", "Importer state count","Exporter log GDP", "Importer log GDP", "log Distance", 
                                               "log Arable Land ratio", "log Population ratio", "log GDPpc ratio"),
                          title = "Concentration of regional military capabilities and dyadic trade. Model 1 includes region, 
                          dyad, exporter, and importer random intercepts. Model 2 includes region, dyad, exporter-year, and importer-year
                          random intercepts. Model 3 includes exporter region, importer region, dyad, exporter, and importer random 
                          intercepts. Note: in models 1, 2, and 4, the exporter region is also the importer region.",
                          label = "table:coefficientsA3",
                          notes.align = "r",
                          star.cutoffs = c(0.05, 0.01, 0.001))

write(x=coef_tableA3, file=paste0(getwd(), "/tables/Coef-TableA3.tex"))


### Appendix Table 4 - regional restrictions
coef_tableA4 <- stargazer(model1.n4, model1.n4.1, model3.n4, model3.n4.1, model5.n4, model5.n4.1,
                          column.labels=c("log Intra-region trade", "log Inter-region trade", 
                                          "Region trade proportion"),
                          column.separate = c(2, 2, 2),
                          model.names=FALSE, style="aer",
                          ci=FALSE, p.auto=TRUE,
                          dep.var.labels.include=FALSE,
                          dep.var.labels="",
                          table.placement="htbp",
                          font.size="footnotesize",
                          model.numbers=TRUE, #float.env="sidewaystable",
                          notes.label=" ",
                          digits=2, notes="", df=FALSE, omit.stat=c("rsq", "aic", "bic", "ll"),
                          order = c(1, 2, 8, 3, 4, 5, 6, 7),
                          omit = "fyear",
                          covariate.labels = c("Regional power NHHI", "Regional power Polity", "NHHI X Reg. power Polity",
                                               "Region cinc total", "Count of states in region", "log GDP", "log Population", "Polity combined score", "Armed conflict"),
                          title = "Replication of Models 1, 3, and 5 from the main text, excluding Western Europe (odd models) 
                          and excluding Western Europe and North America (even models) ",
                          label = "table:coefficientsA4",
                          notes.align = "r",
                          star.cutoffs = c(0.05, 0.01, 0.001))

write(x=coef_tableA4, file=paste0(getwd(), "/tables/Coef-TableA4.tex"))

### Appendix Table 5 - time period restrictions
coef_tableA5 <- stargazer(model1.60, model1.80, model3.60, model3.80, model5.60, model5.80,
                          column.labels=c("log Intra-region trade", "log Inter-region trade", 
                                          "Region trade proportion"),
                          column.separate = c(2, 2, 2),
                          model.names=FALSE, style="aer",
                          ci=FALSE, p.auto=TRUE,
                          dep.var.labels.include=FALSE,
                          dep.var.labels="",
                          table.placement="htbp",
                          font.size="footnotesize",
                          model.numbers=TRUE, #float.env="sidewaystable",
                          notes.label=" ",
                          digits=2, notes="", df=FALSE, omit.stat=c("rsq", "aic", "bic", "ll"),
                          order = c(1, 2, 8, 3, 4, 5, 6, 7),
                          omit = "fyear",
                          covariate.labels = c("Regional power NHHI", "Regional power Polity", "NHHI X Reg. power Polity",
                                               "Region cinc total", "Count of states in region", "log GDP", "log Population", "Polity combined score", "Armed conflict"),
                          title = "Replication of Models 1, 3, and 5 from the main text, examining the Cold War (1960-1988, odd models) 
                          and the post-Cold War period (1989-present, even models)",
                          label = "table:coefficientsA5",
                          notes.align = "r",
                          star.cutoffs = c(0.05, 0.01, 0.001))

write(x=coef_tableA5, file=paste0(getwd(), "/tables/Coef-TableA5.tex"))

### Summary stats
### Note: this table needs some touching up to match the version in the appendix
stat_tableA6 <- stargazer(final.data)


###############
### Figures ###
###############


# Figure 1 - from model 1
m1eff <- effect("NHHI*hegpolity", model1.man, xlevels = list(NHHI = c(0, .2, .4, .6, .8), hegpolity = c(-10, 0, 10)), multiline = TRUE)
forfig1 <- as.data.frame(m1eff)
forfig1$hegpolity.fig[forfig1$hegpolity == -10]  <- "Regional Power Polity = -10"
forfig1$hegpolity.fig[forfig1$hegpolity == 0]  <- "Regional Power Polity = 0"
forfig1$hegpolity.fig[forfig1$hegpolity == 10]  <- "Regional Power Polity = 10"

ggplot(forfig1, aes(x = NHHI, y = fit)) + 
  geom_ribbon(aes(ymin = lower, ymax = upper, linetype=NA), alpha = .8)  + 
  geom_line(size = .5, linetype = 1) + xlab("Regional CINC NHHI") + ylab("log Trade t+1") + 
  # ggtitle("Predicted value of intra-region trade with 95% confidence bounds") +
  facet_wrap(~ hegpolity.fig) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                          linetype="dashed"),
        axis.text.x = element_text(hjust = 1),
        panel.grid.major.x = element_blank(),               
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size=rel(1.8), hjust=0))
ggsave(file=paste0(getwd(), "/figures/figure1.pdf"), width=11, height=10)


# Figure 2 - from model 3
m3eff <- effect("NHHI*hegpolity", model3.man, xlevels = list(NHHI = c(0, .2, .4, .6, .8), hegpolity = c(-10, 0, 10)), multiline = TRUE)
forfig2 <- as.data.frame(m3eff)
forfig2$hegpolity.fig[forfig2$hegpolity == -10]  <- "Regional Power Polity = -10"
forfig2$hegpolity.fig[forfig2$hegpolity == 0]  <- "Regional Power Polity = 0"
forfig2$hegpolity.fig[forfig2$hegpolity == 10]  <- "Regional Power Polity = 10"


ggplot(forfig2, aes(x = NHHI, y = fit)) + 
  geom_ribbon(aes(ymin = lower, ymax = upper, linetype=NA), alpha = .8)  + 
  geom_line(size = .5, linetype = 1) + xlab("Regional CINC NHHI") + ylab("log Trade t+1") + 
  # ggtitle("Predicted value of inter-region trade with 95% confidence bounds") +
  facet_wrap(~ hegpolity.fig) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                          linetype="dashed"),
        axis.text.x = element_text(hjust = 1),
        panel.grid.major.x = element_blank(),               
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size=rel(1.8), hjust=0))
ggsave(file=paste0(getwd(), "/figures/figure2.pdf"), width=11, height=10)


# Figure 3 - from model 5
m5eff <- effect("NHHI*hegpolity", model5.man, xlevels = list(NHHI = c(0, .2, .4, .6, .8), hegpolity = c(-10, 0, 10)), multiline = TRUE)
forfig3 <- as.data.frame(m5eff)
forfig3$hegpolity.fig[forfig3$hegpolity == -10]  <- "Regional Power Polity = -10"
forfig3$hegpolity.fig[forfig3$hegpolity == 0]  <- "Regional Power Polity = 0"
forfig3$hegpolity.fig[forfig3$hegpolity == 10]  <- "Regional Power Polity = 10"


ggplot(forfig3, aes(x = NHHI, y = fit)) + 
  geom_ribbon(aes(ymin = lower, ymax = upper, linetype=NA), alpha = .8)  + 
  geom_line(size = .5, linetype = 1) + xlab("Regional CINC NHHI") + ylab("log Trade t+1") + 
  # ggtitle("Predicted value of regional trade proportion with 95% confidence bounds") +
  facet_wrap(~ hegpolity.fig) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                          linetype="dashed"),
        axis.text.x = element_text(hjust = 1),
        panel.grid.major.x = element_blank(),               
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size=rel(1.8), hjust=0))
ggsave(file=paste0(getwd(), "/figures/figure3.pdf"), width=11, height=10)


## Appendix figures ## 

# Appendix figure 1
m1reff <- effect("NHHI*hegpolity", model1.re, xlevels = list(NHHI = c(0, .2, .4, .6, .8), hegpolity = c(-10, 0, 10)), multiline = TRUE)
forfigA1 <- as.data.frame(m1reff)
forfigA1$hegpolity.fig[forfigA1$hegpolity == -10]  <- "Regional Power Polity = -10"
forfigA1$hegpolity.fig[forfigA1$hegpolity == 0]  <- "Regional Power Polity = 0"
forfigA1$hegpolity.fig[forfigA1$hegpolity == 10]  <- "Regional Power Polity = 10"

ggplot(forfigA1, aes(x = NHHI, y = fit)) + 
  geom_ribbon(aes(ymin = lower, ymax = upper, linetype=NA), alpha = .8)  + 
  geom_line(size = .5, linetype = 1) + xlab("Regional CINC NHHI") + ylab("log Trade t+1") + 
  # ggtitle("Predicted value of intra-trade with 95% confidence bounds") +
  facet_wrap(~ hegpolity.fig) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                          linetype="dashed"),
        axis.text.x = element_text(hjust = 1),
        panel.grid.major.x = element_blank(),               
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size=rel(1.8), hjust=0))
ggsave(file=paste0(getwd(), "/figures/figureA1.pdf"), width=11, height=10)


# Appendix figure 2
m3reff <- effect("NHHI*hegpolity", model3.re, xlevels = list(NHHI = c(0, .2, .4, .6, .8), hegpolity = c(-10, 0, 10)), multiline = TRUE)
forfigA2 <- as.data.frame(m3reff)
forfigA2$hegpolity.fig[forfigA2$hegpolity == -10]  <- "Regional Power Polity = -10"
forfigA2$hegpolity.fig[forfigA2$hegpolity == 0]  <- "Regional Power Polity = 0"
forfigA2$hegpolity.fig[forfigA2$hegpolity == 10]  <- "Regional Power Polity = 10"

ggplot(forfigA2, aes(x = NHHI, y = fit)) + 
  geom_ribbon(aes(ymin = lower, ymax = upper, linetype=NA), alpha = .8)  + 
  geom_line(size = .5, linetype = 1) + xlab("Regional CINC NHHI") + ylab("log Trade t+1") + 
  # ggtitle("Predicted value of inter-region trade with 95% confidence bounds") +
  facet_wrap(~ hegpolity.fig) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                          linetype="dashed"),
        axis.text.x = element_text(hjust = 1),
        panel.grid.major.x = element_blank(),               
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size=rel(1.8), hjust=0))
ggsave(file=paste0(getwd(), "/figures/figureA2.pdf"), width=11, height=10)


# Appendix figure 3
m5reff <- effect("NHHI*hegpolity", model5.re, xlevels = list(NHHI = c(0, .2, .4, .6, .8), hegpolity = c(-10, 0, 10)), multiline = TRUE)
forfigA3 <- as.data.frame(m5reff)
forfigA3$hegpolity.fig[forfigA3$hegpolity == -10]  <- "Regional Power Polity = -10"
forfigA3$hegpolity.fig[forfigA3$hegpolity == 0]  <- "Regional Power Polity = 0"
forfigA3$hegpolity.fig[forfigA3$hegpolity == 10]  <- "Regional Power Polity = 10"

ggplot(forfigA3, aes(x = NHHI, y = fit)) + 
  geom_ribbon(aes(ymin = lower, ymax = upper, linetype=NA), alpha = .8)  + 
  geom_line(size = .5, linetype = 1) + xlab("Regional CINC NHHI") + ylab("log Trade t+1") + 
  # ggtitle("Predicted value of region trade proportion with 95% confidence bounds") +
  facet_wrap(~ hegpolity.fig) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                          linetype="dashed"),
        axis.text.x = element_text(hjust = 1),
        panel.grid.major.x = element_blank(),               
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size=rel(1.8), hjust=0))
ggsave(file=paste0(getwd(), "/figures/figureA3.pdf"), width=11, height=10)


## Figures demonstrating the conditional impact of Polity ##

# Appendix figure 4
m4eff <- effect("NHHI*hegpolity", model1.man, xlevels = list(hegpolity = c(-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10), 
                                                             NHHI = c(.2, .5, .8)), multiline = TRUE)
forfig4p <- as.data.frame(m4eff)
forfig4p$NHHI.fig[forfig4p$NHHI == .2]  <- "Regional CINC NHHI = .2"
forfig4p$NHHI.fig[forfig4p$NHHI == .5]  <- "Regional CINC NHHI = .5"
forfig4p$NHHI.fig[forfig4p$NHHI == .8]  <- "Regional CINC NHHI = .8"

ggplot(forfig4p, aes(x = hegpolity, y = fit)) + 
  geom_ribbon(aes(ymin = lower, ymax = upper, linetype=NA), alpha = .8)  + 
  geom_line(size = .5, linetype = 1) + xlab("Regional power Polity") + ylab("log Trade t+1") + 
  # ggtitle("Predicted value of intra-region trade with 95% confidence bounds") +
  facet_wrap(~ NHHI.fig) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                          linetype="dashed"),
        axis.text.x = element_text(hjust = 1),
        panel.grid.major.x = element_blank(),               
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size=rel(1.8), hjust=0))
ggsave(file=paste0(getwd(), "/figures/figureA4.pdf"), width=11, height=10)

# Appendix Figure 5
m5eff <- effect("NHHI*hegpolity", model3.man, xlevels = list(hegpolity = c(-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10), 
                                                             NHHI = c(.2, .5, .8)), multiline = TRUE)
forfig5p <- as.data.frame(m5eff)
forfig5p$NHHI.fig[forfig5p$NHHI == .2]  <- "Regional CINC NHHI = .2"
forfig5p$NHHI.fig[forfig5p$NHHI == .5]  <- "Regional CINC NHHI = .5"
forfig5p$NHHI.fig[forfig5p$NHHI == .8]  <- "Regional CINC NHHI = .8"

ggplot(forfig5p, aes(x = hegpolity, y = fit)) + 
  geom_ribbon(aes(ymin = lower, ymax = upper, linetype=NA), alpha = .8)  + 
  geom_line(size = .5, linetype = 1) + xlab("Regional power Polity") + ylab("log Trade t+1") + 
  # ggtitle("Predicted value of intra-region trade with 95% confidence bounds") +
  facet_wrap(~ NHHI.fig) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                          linetype="dashed"),
        axis.text.x = element_text(hjust = 1),
        panel.grid.major.x = element_blank(),               
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size=rel(1.8), hjust=0))
ggsave(file=paste0(getwd(), "/figures/figureA5.pdf"), width=11, height=10)


# Appendix Figure 6
m6eff <- effect("NHHI*hegpolity", model5.man, xlevels = list(hegpolity = c(-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10), 
                                                             NHHI = c(.2, .5, .8)), multiline = TRUE)
forfig6p <- as.data.frame(m6eff)
forfig6p$NHHI.fig[forfig6p$NHHI == .2]  <- "Regional CINC NHHI = .2"
forfig6p$NHHI.fig[forfig6p$NHHI == .5]  <- "Regional CINC NHHI = .5"
forfig6p$NHHI.fig[forfig6p$NHHI == .8]  <- "Regional CINC NHHI = .8"

ggplot(forfig6p, aes(x = hegpolity, y = fit)) + 
  geom_ribbon(aes(ymin = lower, ymax = upper, linetype=NA), alpha = .8)  + 
  geom_line(size = .5, linetype = 1) + xlab("Regional power Polity") + ylab("log Trade t+1") + 
  # ggtitle("Predicted value of intra-region trade with 95% confidence bounds") +
  facet_wrap(~ NHHI.fig) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                          linetype="dashed"),
        axis.text.x = element_text(hjust = 1),
        panel.grid.major.x = element_blank(),               
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size=rel(1.8), hjust=0))
ggsave(file=paste0(getwd(), "/figures/figureA6.pdf"), width=11, height=10)




## Dyad model graphs ##

# Appendix figure 7 - Graph for dyad intra-region trade
m1deffr <- effect("NHHI1*hegpolity1", model1dr, xlevels = list(NHHI1 = c(0, .2, .4, .6, .8), hegpolity1 = c(-10, 0, 10)), multiline = TRUE)
forfigA7 <- as.data.frame(m1deffr)
forfigA7$hegpolity1.fig[forfigA7$hegpolity1 == -10]  <- "Regional Power Polity = -10"
forfigA7$hegpolity1.fig[forfigA7$hegpolity1 == 0]  <- "Regional Power Polity = 0"
forfigA7$hegpolity1.fig[forfigA7$hegpolity1 == 10]  <- "Regional Power Polity = 10"

ggplot(forfigA7, aes(x = NHHI1, y = fit)) + 
  geom_ribbon(aes(ymin = lower, ymax = upper, linetype=NA), alpha = .8)  + 
  geom_line(size = .5, linetype = 1) + xlab("Regional CINC NHHI") + ylab("log Exports t+1") + 
  # ggtitle("Predicted value of intra-region dyadic trade with 95% confidence bounds") +
  facet_wrap(~ hegpolity1.fig) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                          linetype="dashed"),
        axis.text.x = element_text(hjust = 1),
        panel.grid.major.x = element_blank(),               
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size=rel(1.8), hjust=0))
ggsave(file=paste0(getwd(), "/figures/figureA7.pdf"), width=11, height=10)



# Appendix figure 8 - two parts
m2deffr <- effect("NHHI1*hegpolity1", model2dr, xlevels = list(NHHI1 = c(0, .2, .4, .6, .8), hegpolity1 = c(-10, 0, 10)), multiline = TRUE)
forfigA8 <- as.data.frame(m2deffr)
forfigA8$hegpolity1.fig[forfigA8$hegpolity1 == -10]  <- "Regional Power Polity = -10"
forfigA8$hegpolity1.fig[forfigA8$hegpolity1 == 0]  <- "Regional Power Polity = 0"
forfigA8$hegpolity1.fig[forfigA8$hegpolity1 == 10]  <- "Regional Power Polity = 10"

f8a <- ggplot(forfigA8, aes(x = NHHI1, y = fit)) + 
  geom_ribbon(aes(ymin = lower, ymax = upper, linetype=NA), alpha = .8)  + 
  geom_line(size = .5, linetype = 1) + xlab("Regional CINC NHHI") + ylab("log Exports t+1") + 
  ggtitle("Exporter region") +
  facet_wrap(~ hegpolity1.fig) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                          linetype="dashed"),
        axis.text.x = element_text(hjust = 1),
        panel.grid.major.x = element_blank(),               
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size=rel(1.8), hjust=0))

m2deffr2 <- effect("NHHI2*hegpolity2", model2dr, xlevels = list(NHHI2 = c(0, .2, .4, .6, .8), hegpolity2 = c(-10, 0, 10)), multiline = TRUE)
forfigA8B <- as.data.frame(m2deffr2)
forfigA8B$hegpolity2.fig[forfigA8B$hegpolity2 == -10]  <- "Regional Power Polity = -10"
forfigA8B$hegpolity2.fig[forfigA8B$hegpolity2 == 0]  <- "Regional Power Polity = 0"
forfigA8B$hegpolity2.fig[forfigA8B$hegpolity2 == 10]  <- "Regional Power Polity = 10"

f8b <- ggplot(forfigA8B, aes(x = NHHI2, y = fit)) + 
  geom_ribbon(aes(ymin = lower, ymax = upper, linetype=NA), alpha = .8)  + 
  geom_line(size = .5, linetype = 1) + xlab("Regional CINC NHHI") + ylab("log Exports t+1") + 
  ggtitle("Importer region") +
  facet_wrap(~ hegpolity2.fig) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust=0, size=rel(1.1)),
        panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                          linetype="dashed"),
        axis.text.x = element_text(hjust = 1),
        panel.grid.major.x = element_blank(),               
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(size=rel(1.8), hjust=0))

# Combine figures 8A and 8B
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  require(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}


pdf(file = paste0(getwd(), "/figures/figureA8.pdf"), width=11, height=10)
multiplot(f8a, f8b, cols = 1)
dev.off()
